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Abstract 



Systems with the coexistence of different stable attractors are widely exploited in 
Q I systems biology in order to suitably model the differentiating processes arising in 

living cells. In order to describe genetic regulatory networks several deterministic 

models based on systems of nonlinear ordinary differential equations have been 

proposed. 

Few studies have been developed to characterize how either an external input 

or the coupling can drive systems with different coexisting states. For the sake of 
>■ ■ simplicity, in this manuscript we focus on systems belonging to the class of radial 

"^ ■ isochron clocks that exhibits hard excitation, in order to investigate their complex 

ly^ ! dynamics, local and global bifurcations arising in presence of constant external 

CN I inputs. In particular the occurrence of saddle node on limit cycle bifurcations is 

t"^ ■ detected. 

o; 
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^ ■ 1. Introduction 



Being the most diffuse formalism to model dynamical systems in science and 
engineering, ordinary differential equations (DDEs) have been widely used in sys- 
tems biology as well. In particular, nonlinear dynamical systems that can exhibit 
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coexisting stable attractors are considered of universal importance, since they al- 
low physiologists to accurately model cell differentiation processes. One of the 
most interesting cases, that goes by the name of of ha rd excitation, is the con- 



currence between oscil latory and non-oscillatory states (|Liul l2002l : iJovanic et al 
20081 : iMinorskvl 1 19741) . 



It is well known the external environment (such as light and temperature) or 
the substrate synthe sis/injection rate play an important role in the dynamics of 
molecular reactions flBastin & Dochainl . Il990 ). In order to properly model these 
effects, ex ternal forcing terms have to be taken into account in the mathematical 
modeling (ILeloup & Goldbeteii 1200 lb . Moreover, the processes inside the cell 
are localized in different spatial domains, while exchanges of chemical species 
take place in the common extracellular medium. Therefore, it is mo re appropriate 
to co nsider reaction diffusio n or multi-compartmental equations ( Jpvanic et al.|, 
2008h . Actually, few studies (tOoldbeter et al.[l200ll : lGoldbeter & Berridge.. 19971) 
have been developed to characterize how an external input or coupling effects can 
drive systems with different coexisting states. 

Typical systems that exhibit such a dynamical behavior are represented by 
the so-called Cyclic Negative Feedback Systems, that arise i n a variety of math- 
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ematical b io-inpired models, from cellul ar signal pathways (iKholodenk 

Liul I2OO2 I ) to gene regul atory networks (ITyson & Othmeii 1 19781 : Ide Jon 

Elowitz & Leibleii |2000|) . Cyclic Negative Feedback systems ( CNF systems) are 
descr ibed by the following set of nonlinear differential equations (|Arcak & Sontag , 
2008h : 

xi = -giixi)+f„{x„) 

Xk = -gk{xk)+fk-i{xk-]), 2<k<n 

wherex^ only assumes positive values, ^p(-) (p= 1, . . . ,n) and/^(-) {q=l,. . . ,n — 
1) are increasing functions, while /„(■) is a decreasing function. The first-order 
time derivative of ;cyt is denoted by Xk- The variables Xk can represent the concen- 
trations of certain molecules (e.g. mRNA, proteins) in the cell. In addition, the 
function /„(•) suitably reproduces the inhibition of xi by the chemical product x„. 
It is easy to deduce from the previous expression that each variable Xk{t) is acti- 
vated by its pre vious neighbor .yt- 1 (t), except for xi (t) that is repressed by x^ (t) . It 
can be proved (|Arcak & Sontagl |2008|) that CNF systems exhibit a global asymp- 



totic stable equilibrium x* = (xj, . . . ,Ji;*), provided the secant criterion (|ThronL 



199 1|) is satisfied . Furthermore, for these typ e of systems the Poincare-Bendixson 
Theorem holds (|Mallet-Paret & Smithlll990|) . that is the coexistence between sta- 
ble equilibria and periodic orbits is allowed. 



In (|Lanza et all 120091) . we have analyzed the effect of coupling on arrays of 
diffusively coupled third order CNF systems. We have shown that CNF arrays 
with diffusive couplings that are constant and local (i.e. they involve only the two 
nearest neighbors of each CNF system) are potentially equivalent to nonlinear 
networks whose elements are fully connected (i.e. each subsystem is linked to 
all the others). Moreover, we have shown that already in a two-compartment 
version of CNF systems new dynamics, such as global periodic oscillations with 
space- variant amplitude (e.g., discrete breathers -like patterns), can arise due to the 
couplings. 

One of the main drawbacks of CNF systems is that they can be of high order 
an d quite difficult to h andle. Even for the network of third order systems analyzed 



in (|Lanza et al.L 120091) . it is matematically complicated to detect and characterize 
all the new periodic solutions to whom the coupling or an external input give 
rise. Therefore, in order to investigate the emergence of such new dynamics, we 
consider a nonlinear dynamical system that is easier to handle and that for certain 
values of its parameters has qualitatively the same dynamical behavior of a CNF 
system. First of all, we focus on the single system and carry out a complete 
analysis of both the local and global bifurcations arising in presence of constant 
external inputs. In particular, we show the occurrence of saddle node on limit 
cycle bifurcations and that in presence of such bifurcations our system behaves as 
a relaxation oscillator. 

This manuscript is structured as follows. In Section 2 we introduce the model 
under study, i.e. a radial isochron clock with hard excitation that is well known 



in literature as the normal form for the Bautin bifurcation (|Kuznetsovl |2004j; 



Izhikevic h. 2001). Moreover, we investigate how this system changes its form 



due to the presence of a constant external input. In Section 3 a complete analysis 
of local and global bifurcations arising in presence of constant external inputs is 
carried out. In particular, we show that the periodic solutions of our system can 
disappear through saddle node on limit cycle bifurcations. Section 4 is devoted to 
conclusions. 

2. The case of radial isochron clocks 

Let us consider a simpler model with respect to CNF systems but with the 
same dynamics: 

z = {oo + jao)z + oi\z\z + 02\z\^z + lo, (1) 



where z = re,^'^ G C is a complex variable, f2o > and Oq, Oi , 02 are real parame- 
ters. The term /q represents the action of a constant external input and from now 
on we will assume /q > without any loss of generalitjo 

Absence of constant external input 

In the (r,(|)) coordinates, system ([T]) with /q = can be recast as: 



r = r(oo + 0ir + 02r^ 



(2) 



It is worth noting tha t Q belongs to the category of the radial isochron clocks 
proposed by Winfree (|Winfreel . l200l|) . and it h as been widely exploited a s a sim- 
plified model for the Hodgkin-Huxley neuron (IBoushel & CurraiJ, l2007h and for 
studying the phenomenon of cardiac fibrillation ddePaoii 'l994). 

This system can be easily investigated since the study of the periodic solutions 
of © reduces to the analysis of the equilibria of the first equation of ([21). It is easy 
to see that 



a Hopf bifurcation takes place for Oq 
and supercritical otherwise; 



0, which is subcritical for Oi > 



a double limit cycle bifurcation (also known as fold, or tangent, or saddle- 
node bifurcation of limit cycles) occurs for Oj — AoqOi = 0, and Oi > 0. 
A double limit cycle bifurcation occurs when a branch of stable periodic 
solutions and a branch of unstable perio dic solutions coalesce and oblit- 



erate each other at the bifurcation point (|Guckenheimer & Holmesl 11983 
Kuznetsovll2004 : 



this system can undergo a Bautin bifurcation, when a H opf bifurcation and 
a double l i mit cy cle bifurcation occur simultaneously (|Kuznetsovl 12004 : 
Izhikevichl 1200 1|) . In our case it happens for Oo = Oi = but 02 ^ 0. When 



02 < 0, the Bautin bifurcation is said to be supercritical, while for 02 > 
it is subcritical. In particular, for 02 < the cycle with larger amplitude is 
stable. 

The complete bifurcation diagram is represented in Figure [B 



It is easy to notice that the system in invariant under the transformation z — > — z,/o — > —to- 
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Figure 1 : Bifurcation diagram for the single radial isochron clock ^. The different bifurcation 
curves are labeled as follows: Hi stands for subcritical Hopf bifurcation, H2 for supercritical Hopf 
bifurcation, DC for double limit cycle bifurcation, and B for Bautin bifurcation. The dotted line 
where a DC bifurcation occurs has equation a\ — O0O2 = 0. 



In the following, we are interested in studying a system that presents a hard 
excitation behavior, thus we focus on a choice of the parameters such that we have 
the coexistence of a stable equilibrium point and a stable limit cycle: 



(3) 



In Figure [2] the phase portrait of our system is represented. We have a sort of 
concentric structure of cycles, where the origin is a stable equilibrium point, sor- 
rounded by two closed curves alternatively stable and unstable. Thus, it is easy 
to understand wh y a similar configuration goes by the name of hard excitation 
(IMinorskyl 119741) : if the systems is in the steady state, then it needs a strong 
perturbation to cross the separatrix and start to oscillate. 




Presence of external constant input 

In order to simplify the notation and reduce the number of control parameters, 
we rescale the variables involved in ([T]) in the following way: 



U,0^ 




|02| |Oo| 

Thus, from ([T]) we obtain the system 

z = (sgnao + j^)z + a|z|z + sgna2|z|^z + /, 



(4) 




Figure 2: Phase portrait of a system with hard excitation. The stable equilibrium point and the 
stable limit cycle (solid line) are separated by an unstable limit cycle (dashed line). 



where sgn(-) is the sign function and 
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(5) 



Conditions ([3]) imply sgnoo = — l,sgn02 = — l,a > 2, and thus system dH) be- 
comes 

z=i-l+j^)z + a\z\z-\z\^z + I. (6) 

Exploiting the cartesian coordinates (z = x + jy), system © can be rewritten as 



x= {-1+a^/x^ +y^ - {x^ +y^))x-Q.y+I 
y = {-\ +a^/x^ +y^ - {x^ +y^))y + ax, 



or to simplify notation 



x = g{x,y)x-Q.y+I 



(7) 



(8) 



where ^(.x:,j) = {-1 +a^/x^ +y^ - [x^ + y^)). 

First of all, in order to characterize the dynamics of this system, we are in- 
terested in finding all the equilibrium configurations. Thus, we have to solve the 
following set of nonlinear algebraic equations: 



g{x,y)x-Qy + I = 

g{x,y)y + ^ = o. 



(9) 



From the second equation we obtain 






(10) 



having assumecO J 7^ 0- The substitution of expression (flOl) in the first equation 
of ^ yields 



X 



-ri D.y + 1 

y 



0. 



(11) 



The solutions of ([8]) are precisely the intersection points between the circumfer- 
ence ^ ■ x^ + {y— '^) =(2^) ^^'^ ^^^ curve S : X = —j^yg{x,y). Moreover, 
from (fTTI) we can conclude that all the equilibrium configurations have a positive 
J— component. Furthermore, since from (fTTI) we have x^+y^ = -^y, we can notice 
that actually S has the following expression: 



S : X 



1 

a 



l+a 




a^ Q. 



y 



and we can derive that the equilibrium points of system ([8]) are the intersections 
between the curves: 

v2 / I \2 



'^■■-' + iy-^) =i^y 



x= — 



Q. 



'^+^\fky-iy 



(12) 



Substituting the second equation of (fT2)) in the first one, and introducing the new 

variable C, = \ -^y, we finally conclude that the equilibrium points of system dV]) 
are the roots of the following polynomial: 



^^-2d:J' + {o? + 2%'^-2a:,^ + {\+Q}%^-I^ = Q. 



(13) 



It is worth noting that, due to the definition of ^, we are interested only in the real 
and positive roots of (fT3] ). 

By choosing a = 3, in Figure |3] the equilibrium configurations and their sta- 
bility properties as function of / for different fixed values of Q are represented. 
Depending on the values of the two parameters ^ and /, a different number of 
solutions and therefore different dynamical behaviors are admissible. 



^The case y = Q can be treated separately and is not interesting for our study. In fact, from (|9]l 
we can conclude that y = implies x = and 7 = 0. Thus, a solution with y = Q can be achieved 
only in absence of external input. 
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Figure 3: Stable (in blue) and unstable (in red) equilibrium configurations of system ^ as function 
of the parameter I, for different values of ii (in (a) and (b) i2 = 0.25, in (c) and (d) ii = 0.75, in 
(e) and (f) H = 1 .5). Here a is set equal to 3. In the left and right columns the x and y-components 
are represented, respectively. The bifurcations are indicated as follows: SN= saddle-node, SNLC= 
saddle-node on limit cycle, H= Hopf. 



3. Bifurcation analysis 

In this section we study the bifurcations occurring in system (|7]) as the param- 
eters {0.,I) are varied. 

We have seen in the previous section that in absence of external input (/ = 0) 
the system presents two nested periodic solutions, one stable and one unstable, 
for every value of Q.. For small values of /, these solutions are mantained, but, 
increasing the intensity of the input, they disappear through a sequence of different 
both local and global bifurcations. 

3.1. Local bifurcations 

In order to investigate local bifurcations of equilibria in system (|7]) we shall 
look at the polynomial discriminant of (fT3l) . and at the linearization of (|7]) in the 
neighborhood of the equilibrium points. A polynomial discriminant is defined as 
the product of the squares of the differences of the polynomial roots Sj. For a 
polynomial of degree n in the form 

q{x) = a„x" +an-\:x"^ H \-aix + ao — 



the discriminant is defined as (|Cohenl . Il993|) 



n 



Since the discriminant vanishes in presence of a multiple root, the values of Q. 
and / for which the discriminant D is equal to zero identify loci of coalescences 
or births of solutions, and therefore possible bifurcations. In Figure |4] the curves 
Z) = for (fT3l) in the Q. — I plane for a = 3 are represented. 

Now we turn our attention at the linearization of ([7]). The Jacobian matrix of 
system Q, evaluated in the equilibrium points {x,y), has the following expression: 



g{x,y) +x^h{x,y) xyh{x,y) - Q 
xyh{xj)+Q. g{xj)+y^h{x,y) 



J _ I SK-^iJJ ^-^ '^K-^:j; ^y-K-^TJI ^' \ /^c\ 



where 



h{x,y) = —^=-2. (16) 
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Figure 4: The locus of parameters where the discriminant of polynomial (fTsT i with a = 3 annihi- 
lates. This points out the possible occurrence of bifurcations of equilibria in the system. 

The sign of the real part of the eigenvalues can be determined by looking at the 
trace and the determinant of the Jacobian matrix. In the present case they are given 
by 



_2 , -2\i(- -\ 



ixJ = 2g{x,y) + {x^ + y^)h{xj) 



^2 ltt2 tjl L „2^^ ^(^2 ,-2\ 



- l{-l+asJT+y^ - {x^+y^)) + {T +y 



2(-l+aC-CVocC-2C^ 

.^2 



a 



^x^+i^ 



= -4C +3a^-2, 



(17) 



diaJ={g{x,y)+x^h{x,y)){g{x,y)+y^h{x,y))-{xyh{x,y)-^){xyh{x,y) + Q.) 
= {g{x,y)f + {x^ + f)g{x,y)h{x,y)-Q} 

ni _ ,__, 

= -^ — Ixh{x,y) 



ai I 

y £2 



f ^ 



'^''yy-y 



a 



= ^ + C|-l+aC-fl(cx-2g, 



(18) 



having exploited relations (|TOl) . ([T2|) . and (fT6l) . and having introduced the ancillar 



10 



variable C, = 

Time-domain simulations show that for certain values of Q. and / system Qi 
exhibits Hopf bif urcations. The con ditions for the occurrence of this bifurcation 
are the following (|Kuznetsovl |2004|) : 




From (fT71 ). condition tr7 = leads to 



■4^%3aC-2 



0, 



that is 



= 1,2 



3a±V9a2-32 
8 ■ 



(19) 



(20) 



(21) 



These two solutions exist only for 9a^ - 32 > 0, that is only for a > ^ = 1.886. 
However, this condition is alway satisfied, since a > 2 due to ([3]). 

From the definition of C, and from (fT2l) . we can find the coordinates of the equi- 
librium points Pi = (lij^i) and P2 = (12,^2)' ^t which the trace of the Jacobian 
annihilates: 



yi 



Xi 



1 



^ ^2 



l+oc^z-C, 



/=1,2. 



Substituting (|2T)) into (fT3l) we obtain the Hopf bifurcation curves in the parameters 
plane 



■7^ 



slA 



I=lf = ^i -2aC, +(a^ + 2)C, -2aC, +(1+^^^)^- /= 1,2. (22) 
For instance, with a = 3 (the case considered in Figure [3]) we have 



Ci 



Pi 



Q. 



Pi 



16a 16/1 

25 + 256a2 



4096 
1 AQ.\ 



7f = 4(1+^2 



11 



Actually , we still have to discriminate between Hopf bifurcations and neu- 
tral saddles (|Kuznetsovl |2004|) . exploiting the condition on the determinant of the 

Jacobian matrix. Introducing again the variable C,i = y^y/, and recalling the ex- 
pression for /, in (|22)) . we get 



det7| 



li 
Requiring the positivity of det7, it yields 



det7|^ > 



^2. .74 , C..73 o/.,2 -^-2 



a' > -3^i +5aC,- -2(a^ + 2)^,. +3a^,- 1 (23) 



64 



-a" 



1 \ 



aK,+ 



4 / 



1 1 

4~32' 



-a 



where we used (|20l ) to reduce the degree of the right hand side. In particular, 
substituting the values of ^, found in (|2TI) . we obtain 



det7|^^ > 



det7|^^ > 



^2 9 4 1 1 2 
Q. > a -I a 

^512 ^4 8 



0-2-' -h^)^^^^' 



(24) 



^2 9 4 1 1 2 / 3 3 



512 



4 8 



512 



^aj V9a2-32. 



Finally, we can conclude that for a = 3 in the plane Q. — Ithe Hopf bifurcations 
occur if the following conditions are satisfied (see Figure [5]): 



25+256^-' 



4096 



Q > — 
li ^ 16 



1^ = 4(1+0^) 

a> 1. 



(25) 



It is worth observing that, for the examples shown in Figure [H we have one 
Hopf bifurcation for Q. = 0.75 (in this case only Pi exists) and two Hopf bifur- 
cations for Q. = 1.5. These situations are in perfect agreement with conditions 
(125]). 
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Figure 5: Curves of Hopf bifurcations in the Q. I plane for a = 3. 

Remark 1. The points Pi and P2 found above have a direct relation with the 
curves T and S in (IT2l) . If we consider S as S : x = f{y), then it is possible to 
show that P\ and P2 are the relative maximum and minimum of S, respectively. 
Furthermore, the values I\ and I2 are precisely the ones for which these relative 
maximum and minimum points lie on the circumference T. 



To determine the stability of the emerging lim it cycle, we use normal form 
theory. We recall (|Guckenheimer & HolmesLll983h that for a planar system in the 
form 

with F{0) = G(0) = and DF{0) = DG{0) = 0, we have to evaluate the following 
quantity: 



_ 1 

1 



(27) 



+ -j-7— {Fxy{Fxx + Fyy) — Gxy{Gxx + G 



yy) 



' fxx^xx ~T PyyLr 



yy^yyj , 



where all the partial derivatives are computed in (0,0). Thus, if a < we can 
conclude that the Hopf bifurcati on is supercritical, while if a > is subcritical 
(|Guckenheimer & Holmeslll983i) . 

In order to proceed with the computation we have to move the bifurcation 
points from Pi and P2 to the origin. Let us introduce the following change of 



13 



variables: 

"=-'~^' (28) 

v = y-yi, 

where / = 1,2, depending on the point we are interested in. Thus, our system ^ 
can be recast as 

where 

Fiu,v)= l~l+a^{u + Xi)^ + {v + yi)^-{{u + Xi)^ + {v + yif)]{u+Xi)-Q.yi+I 

G{u,v)= l-l+a^J{u + Xi)^ + {v + yi)^-{{u + Xi)^ + {v + yi)^)]{v+yi) + Q:xi. 

It is easy to check that the functions F{u,v) and G{u,v) satisfy the conditions 
above. Furthermore, evaluating their partial derivatives, it is possible to observe 
that Gu = Fy, which implies Guu = F^y, Guv = ^w, and Guuv = Fuw Expression 
(t27l) reduces to 

a = -z {Fuuu + ^Fuvv + Gvvv) • (30) 

Id 

Computing these partial derivatives and evaluating them in (0,0), we get 

Recalling that from (fTTl) we have xj+yj = ^y, = ^, , we obtain 

a=^(-16 + |) (32) 

and therefore 

1 / 3a\ 3a 

a, = -16 + — =-1 + >0. (33) 

^'^ 16 V Ciy 2(3a-V9a2-32) 

Therefore, we conclude that for every a > 2 and Q. > j^, in Pi we have a subcrit- 
ical Hopf bifurcation. Analogously, evaluating the quantity a in P2 we obtain: 

1 / 3a\ 3a 
a, =— -16 + ^ =-1 + <0. (34) 

^'^ 16 V ^2 J 2(3a+V9a2-32) 

14 



In this case we expect a supercritical Hopf bifurcation for every a > 2 and ^ > 1. 
Both the resuks have been confirmed by time-domain numerical simulations (see 
also Figure |3]). 

We now turn our attention to saddle-node bifurcations of equilibria, that are 
characterized by d et7 = 0, since they involve the presence of a null eigenvalue 
(|Kuznetsovll2004|) - Using ^T^ we obtain 



/2 = -2^% 3a^^ - (2 -f a2)f + a^l 



(35) 



Because the Jacobian is evaluated in the equilibrium configurations, we recall 
that C, will be a root of (fT3l) . Thus, substituting (l35l) in (fT3l) we obtain 



-7;4 



^ 3^ -5aC+2(2 + a^)^ -3a^+(l+^' 



0. 



(36) 



We conclude that, for any fixed value of the parameter D., we have a saddle-node 
bifurcation at the equilibrium points that satisfy: 



—4 

3C • 



■5aC +2(2 + a2)^ -3aC + (l+fl2 







for the values of 7 given by (I35I) . 

Remark 2. The conditions for the occurrence of a saddle-node bifurcation leads 
to the same set of curves in the Q. — I plane of Figure^ This is due to the fact that 
a zero eigenvalue for the Jacobian matrix of a generic dynamical system x = f{x) 
implies that the equilibrium poin t has a multiplicity equal to two as zero of the 
function f(x)=0 AKuznetsou \2004^) . 



Finally, we consider the case 




(37) 



that corresponds to a codimension 2 bifurcation, the so-called Bogdanov-Takens 
one. In our case, we have previously seen that the trace of the Jacobian matrix is 
equal to zero at the points P\ and P2, with /,- given by (|22|) . The further condition 
on the determinant of 7 leads to the following critical values of Q.: 

9 



Pi 
Pi 



Q.i 



J^? = 



512 

9 

512' 



-a 



-a 






(38) 
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In particular for a = 3 we have 

^^ , 5 a 



Pi 



16^'16/i 



1 4a 



with /i 



V25 + 256^2 
8 



._,_) with /2 = 2Vl+^2 



and, exploting (|38l) . we conclude that in 

{aJ)=(^,^^f2\ and (a,/) = (l,2v^ 
we have two Bogdanov-Takens bifurcations. 



3.2. Global bifurcations 

Time-domain numerical simulations reveal that, for some couple (a,/) shown 
in Figure HI the saddle-node bifurcations of equilibria are nontrivial, in the sense 
that they involve the appearance and disappearance of limit cycles. 

Firstly introduced and studied in (lAndronov et al.L Il973h . in system (|7]) it in- 
volves the disappearance of the periodic solutions. In lite rature, this bifurcation 
goes by several names: saddle-node on limit cy cle (SNLC ) mopp ensteadt & Izhikevich , 
Il997.) . saddle-node on i nvariant cycle (SNIC) dlzhikevich. 2006), saddle-node in- 
finite peri od (SNIPER) (IMcCormick et al.L Il99l|) or saddle-node homoclinic bi- 
furcation (JKuznetsovL l2004h . In particular, its name saddle-node on invariant cycle 
is due to the fact t hat it is a standard saddle-node, but occurs on an invariant cycle 
(jIzhikevichL 12006 1). Let us suppose that the system exhibits a periodic solution, as 
in Figured (a). The emergence of a saddle-node (see Figure |6](b)) coincides with 
the break of the limit cycle, that becomes a homoclinic trajectory. As the bifur- 
cation parameter increases (see Figure |6](c)), the node and the saddle move away 
each other and two heteroclinic trajectories arise to connect the two equilibria. 

It is worth remarking that this bifurcation is both gl obal and local, as i t involves 
a simultaneous collision of equilibria and manifolds (IKuznetsovi |2004() . In fact, 
as we have shown in the previous section, local analysis describes only the saddle- 
node bifurcation of equilibria, missing the disappearance of the limit cycle. 

In general, detecting a homoclinic trajectory is not a simple task. A possible 
evidence of a SNLC bifurcation is given by the period of the limit cycle. In fact, 
the closer is the parameter to the critical value, the larger is the period of the cor- 
responding limit cycle, that tends to infinity approaching the SNLC bifurcation. 
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(a) (b) (c) 

Figure 6: Saddle-node bifurcation on limit cycle (SNLC). 
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Figure 7: Period T of the unstable (a) and stable (b) limit cycles as function of the external input 
/, for Q. = 0.25. In Id = 0.091 and Id — 0.653 the system undergoes two SNLC bifurcations, that 
entail the disappearance in succession of the two limit cycles. 



Furthermore, the system spends more time near the place where the saddle-node 
will appear, in a sort of sense having a hunch of the future bifurcation. 

As an example, let us consider the case Q. = 0.25, for which it is possible to 
see that our system undergoes two SNLC bifurcations. In Figure |7] the periods 
of the two (stable and unstable) limit cycles as function of the external input are 
shown. It is worth observing that the two periods tend to infinity, approaching the 
bifurcation points Id = 0.091 and /c2 = 0.654 (see also Figures Ha-b)), respec- 
tively. 

In Figure [8] the wave forms for the x-component of the two limit cycles in 
proximity of the respective SNLC bifurcations are represented. Notice the analogy 
with dynamics of the so-called relaxation oscillators. 

It is possible to see that the system tends to stay for the majority of the time 
around one configuration, that is precisely where the saddle-node appears at the 
bifurcation point. In order to characterize the positions in the phase plane where 
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8: Wave form of the unstable (a) and stable (b) limit cycles (x-component), in proximity of 
bifurcations (/ = 0.090 and / = 0.654, respectively). The frequency Q. is set equal to 1. 



the system spends most of the time, a histogram of the dynamics arising from the 
simulations in the time-domain {tfi„ = 1500 and At = 0.01) has been computed 
(see Figure (9]). Since for Id = 0.091 the saddle-node has coordinates SNl ~ 
(0.1081,0.3284), while for Ic2 = 0.654 we have SNl = (-0.1368,2.6066), the 
check with Figure [9] leads to the expected conclusions. The same approach has 
been carried out for both the two limit cycles in absence of external input (/ = 0), 
to make a comparison with the previous relaxation oscillator-like behavior. In this 
case, basically all the states are uniformly visitated by the sistem. The two peaks 
at minimum and maximum values are due to the discretization of the variables for 
the histogram computation. 

It is worth observing that, in order to obtain Figures |71 [8] and [9] for the unstable 
limit cycle, the numerical simulations have been performed back in time. 

4. Conclusion 

Cyclic Negative Feedback Systems (CNF systems) are one of the most ex- 
ploited mathematical frameworks to model phenomena arising in systems biology, 
such as cascades of molecular reactions inside the cell. Since it is well known that 
these events take place in different compartments, it seems more appropriate to 
consider networks of diffusively coupled CNF systems. In addition, the effect of 
an external agent, such as the light, the temperature or the synthesis of the sub- 
strate, is suitably modeled as a constant external term. Unfortunately, due to the 
huge number of factors involved in such mechanisms, networks of coupled CNFs 
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Figure 9: Histogram of the unstable (a-c) and stable (b-d) limit cycle dynamics near the SNLC 
bifurcation (/ ~ 0.090 and / = 0.653 in (a) and (b), respectively) and far away from the critical 
values of the bifurcation parameter (/ = in both (c) and (d)). 



systems can be high-dimensional, and therefore their dynamics can be difficult to 
fully characterize. 

Inspired by the complex spatio-temporal patterns found in arrays of diffusively 
coupled Cyclic Negative Feedback systems (CNF systems), we have considered 
the case of radial isochron clocks that exhibit the coexistence of different stable 
attractors, as well as CNF systems. In fact, these systems present the same qual- 
itative behavior of CNF systems, but they are quite easy to handle since they are 
of lower order. In particular, we have dealt with systems that present a hard exci- 
tation behavior, i.e., that display at the same time a stable equilibrium point and a 
stable limit cycle. 

As a first step, we have focused on the effect of a constant external input on 
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a single radial isochron clock with hard excitation. We have carried out a charac- 
terization of local and global bifurcations, and in particular, we have detected the 
occurrence of saddle-node on limit cycle bifurcations. It is interesting to notice 
that, in presence of such bifurcations, the system exhibits a relaxation oscillator- 
like dynamics. 

Once we have completely analyzed the dynamical behavior of a single system 
in presence of a constant external input, we are now interested in considering 
networks of diffusively coupled radial isochron clocks, in order to investigate the 
complex dynamics that may arise due to the additional effect of coupling. 
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